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Abstract. Numerical analysis has no satisfactory method for the more realistic 
optimization models. However, with constraint programming one can compute 
a cover for the solution set to arbitrarily close approximation. Because the use 
of constraint propagation for composite arithmetic expressions is computation- 
ally expensive, consistency is computed with interval arithmetic. In this paper we 
present theorems that support, selective initialization, a simple modification of 
constraint propagation that allows composite arithmetic expressions to be han- 
dled efficiently. 

1 Introduction 

The following attributes all make an optimization problem more difficult: having an 
objective function with an unknown and possibly large number of local minima, being 
constrained, having nonlinear constraints, having inequality constraints, having both 
discrete and continuous variables. Unfortunately, faithfully modeling an application 
tends to introduce many of these attributes. As a result, optimization problems are 
usually linearized, discretized, relaxed, or otherwise modified to make them feasible 
according to conventional methods. 

One of the most exciting prospects of constraint programming is that such difficult 
optimization problems can be solved without these possibly invalidating modifications. 
Moreover, constraint programming solutions are of known quality: they yield inter- 
vals guaranteed to contain all solutions. Equally important, constraint programming 
can prove the absence of solutions. 

In this paper we only consider the core of the constraint programming approach to 
optimization, which is to solve a system of nonlinear inequalities: 

gi(xi , x% , . . . , x m ) < 

g 2 {xi ,x 2 ,..., x m ) < ^ 

9k{xi , X2 , . .. , x m ) < 

It is understood that it may happen that g; = — gj for some pairs i and j, so that equali- 
ties are a special case. If this occurs, then certain obvious optimizations are possible in 
the methods described here. 
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The ability to solve systems such as supports optimization in more ways than 
one. In the first place, these systems occur as conditions in some constrained optimized 
problems. Moreover, one of g±, . . . , could be defined as f{x\ , . . . , x m ) — c, where / 
is the objective function and where c is a constant. By repeatedly solving such a system 
for suitably chosen c, one can find the greatest value of c for which Q is found to have 
no solution. That value is a lower bound for the global minimum 1151 . 

This approach handles nonlinear inequalities with real variables. It also allows some 
or all variables to be integer by regarding integrality as a constraint on a real variable 


All constraint programming work in this direction has been based on interval arith- 
metic. The earliest work |6| used a generic propagation algorithm based directly on 
domain reduction operators for primitive arithmetic constraints. These constraints in- 
cluded sum(x, y, z) defined as x + y — z for all reals x, y, and z. Also included was 
prod(x, y, z) defined as xy — z for all reals x, y, and z. 

This was criticized in 1 3 1 which advocated the use of composite arithmetic expres- 
sion directly rather than reducing them to primitive arithmetic constraints. In |4 14 1 it 
was acknowledged that the generic propagation algorithm is not satisfactory for CSPs 
that derive from composite arithmetic expressions. These papers describe propagation 
algorithms that exploit the structure of such expressions and thereby improve on what is 
attainable by evaluating such expressions in interval arithmetic. Selective Initialization 
was first described in 1 16 1. This was done under the tacit assumption that all default do- 
mains are [—00, +00]. As a result some of the theorems in that paper are not as widely 
applicable as claimed. 

All of these researches are motivated by the severe difficulties experienced by con- 
ventional numerical analysis to solve practical optimization problems. They can be re- 
garded as attempts to fully exploit the potential of interval arithmetic. In this paper we 
also take this point of view. We show that, though Equation Q can contain arbitrar- 
ily large expressions, only a small modification of the generic propagation algorithm is 
needed to optimally exploit the structure of these expressions. This is made possible by 
a new canonical form for Equation that we introduce in this paper. In addition to sup- 
porting our application of constraint processing to solving systems similar to Equation 
0, this canonical form exploits the potential for parallelism in such systems. 



2 A software architecture for optimization problems 

In Figure^we propose a hierarchical software architecture for optimization problems. 
Each layer is implemented in terms of the layer below. 

In the introduction we briefly remarked on how layer 4 can be reduced to layer 3. 
More detail is given in 1 15 1. For the transition between layers and 1 there is much ma- 
terial in the interval arithmetic literature. The part that is relevant to constraint process- 
ing can be found in 1 9 1 . In the present paper we present a new method for implementing 
layer 3 in terms of layer 2. But first we review the transition between layers 1 and 2. 
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Floating-point arithmetic 
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Hardware 



Fig. 1. A software architecture for optimization problems. 

3 Preliminaries 

In this section we provide background by reviewing some basic concepts. These reviews 
also serve to establish the terminology and notation used in this paper. The first few 
sections apply to all constraint satisfaction problems, not only to numerical ones. 

3.1 Constraint satisfaction problems 

A constraint satisfaction problem (CSP) consists of a set of constraints. Each of the 
variables in the constraint is associated with a domain, which is the set of values that 
are possible for the variable concerned. Typically, not all sets of values can be domains. 
For example, sets of real values are restricted to intervals, as described later. 

A valuation is a tuple indexed by variables where the component indexed by v is an 
element of the domain of v. A solution is a valuation such that each constraint is true if 
every variable in the constraint is substituted by the component of the valuation indexed 
by the variable. The set of solutions is a set of valuations; hence a set of tuples; hence 
a relation. We regard this relation as the relation defined by the CSP. In this way the 
relation that is the meaning of a constraint in one CSP can be defined by another. This 
gives CSPs an hierarchical structure. 

With each constraint, there is an associated domain reduction operator, DRO for 
short. This operator may remove from the domains of each of the variables in the con- 
straint certain values that do not satisfy the constraint, given that the other variables 
of the constraint are restricted to their associated domains. Any DRO is contracting, 
monotonic, and idempotent. 

When the DROs of the constraints are applied in a "fair" order, the domains con- 
verge to a limit or one of the domains becomes empty. A sequence of DROs activations 
is fair if every one of them occurs an infinite number of times. The resulting Cartesian 
product of the domains becomes the greatest common fixpoint of the DROs II 3111 . If 
one of the domains becomes empty, it follows that no solutions exist within the initial 
domains. 
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In practice, we are restricted to the domains that are representable in a computer. As 
there are only a finite number of these, any fair sequence of DRO applications yields 
domains that remain constant from a certain point onwards. 

3.2 Constraints 

According to the usual terminology in constraint programming, a constraint states that a 
certain relation holds between its arguments. But in first-order predicate logic the same 
role is played by an atomic formula. In this paper we adopt the terminology of first- 
order predicate logic for the meaning of "atomic formula" and we reserve "constraint" 
for a special case. 

Thus an atomic formula consists of a predicate symbol with terms as arguments. 
A term is a function symbol with terms as arguments. What makes an atomic formula 
first-order is that a predicate symbol can only occur as the outermost symbol in the 
formula. 

At first sight, the inequalities in Equation Q are atomic formulas. In fact, they 
follow the related, but different, usage that prevails in informal mathematics. The in- 
equality 

gi(xi , x 2 , ■ ■ ■ , x m ) < (2) 

is an atomic formula where < is the predicate symbol with two arguments, which are the 
terms gi(x\, X2, ■ ■ ■ , x m ) and 0. A possible source of confusion is that in mathematics 
9i(x%, X2, ■ ■ ■ , x m ) is not necessarily interpreted as a syntactically concrete term, but 
as an abstractly conceived function g± applied to the arguments x±, . . . , x m . The func- 
tion could be defined by means of a term that looks quite different; such a term could 
be nested and contain several function symbols. For example, Equation l|2) could be 
g%{x\, X2) < with gi defined as gi(x, y) — x 2 + xy — y 2 for all x and y. Accordingly, 
the atomic formula corresponding to Equation (|2jl is 

< (+(sq(x),-(x(x,y),sq(y))),0). (3) 

Taking advantage of infix and postfix notation this becomes x 2 + xy — y 2 < 0. 

A constraint is an atomic formula without function symbols. An example of such an 
atomic formula is sum(x, y, z), which is a ternary constraint whose relation is defined 
by x + y = z for all reals x, y, and z. 

In this paper we translate Equation to a CSP with the set of constraints 

{ti < 0, sum(t 2 ,t 3 ,ti), sq(x, t 2 ), sum(t 5 , t 3 , U),prod(x, y, t s ), sq{y, i 4 )} 

Consider a constraint c(x±, . . . , x n ). The meaning of predicate symbol c is a rela- 
tion, say r. For all i 6 {1, . . . . n}, a value a, for variable Xi is inconsistent with respect 
to r and domains X\, . . . , Xi + i, . . . , X n iff it is not the case that 

3ai e Xi, . . . , 3a 4 _i e JQ-i, 3a i+ i e X i+ i, . . . , 3a n G X n such that 

(ax,...,a n ) e r. 
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A DRO for c may replace one or more of the domains Xi, . . . , X n by a subset of it 
if the set difference between the old and the new domains contains inconsistent values 
only. A DRO is optimal if every domain is replaced by the smallest domain containing 
all its consistent values. We call a constraint primitive if an optimal DRO is available 
for it that is sufficiently efficiently computed. What is sufficient depends on the context. 

3.3 Constraint propagation 

To gain information about the solution set, inconsistent values are removed as much as 
possible with modest computation effort. For example, DROs can be applied as long 
as they remove inconsistent values. It is the task of a constraint propagation algorithm 
to reach as quickly as possible a set of domains that cannot be reduced by any DRO. 
Many versions of this algorithm exist 1 1/211. They can be regarded as refinements of the 
algorithm in Figure|2] which we refer to as the generic propagation algorithm; GPA for 
short. 

GPA maintains a pool of DROs, called active set. No order is specified for applying 
these operators. Notice that the active set A is initialized to contain all constraints. 

put all constraints into the active set A 
while ( A ^ 0) { 

choose a constraint C from A 

apply the DRO associated with C 

if one of the domains has become empty, then stop 

add to A all constraints involving variables whose domains have changed, if any 
remove C from A 

} 

Fig. 2. A pseudo-code for GPA. 



3.4 Intervals 

A floating-point number is any element of F U {— oo, +oo}, where F is a finite set 
of reals. A floating-point interval is a closed connected set of reals, where the bounds, 
in so far as they exist, are floating-point numbers. When we write "interval" without 
qualification, we mean a floating-point interval. For every bounded non-empty interval 
X, lb(X) and rb(X) denote the least and the greatest element of X respectively. They 
are referred to as the left and the right bound of X. If X is not bounded from below, then 
lb{X) = — oo. Similarly, if X is not bounded from above, then rb(X) = +oo. Thus, 
X = [lb(X),rb(X)] is a convenient notation for all non-empty intervals, bounded or 
not. 

A box is a Cartesian product of intervals. 
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3.5 Solving inequalities in interval arithmetic 

Moore's idea of solving inequalities such as those in Equation Q by means of interval 
arithmetic is at least as important as the subsequent applications of interval constraints 
to this problem. 

Suppose we wish to investigate the presence of solutions of a single inequality in 
Equation in a box X\ x • • • x X m . Then one evaluates in interval arithmetic the 
expression in the left-hand side. As values for the variables x±, . . . , x m one uses the 
intervals X\, ■ ■ • , X m . Suppose the result is the interval [a^, bi]. We have exactly one 
of the following three cases. If < a, for at least one i, then there are no solutions. If 
bi < for all i, then all tuples in Xi x • • • x X m are solutions. If a,; < < bi for at 
least one i, then either of the above may be true. In this same case of a; < < bi, it 
may also be that some of the tuples in X\ x • • • x X m are solutions, while others are 
not. 

Again, in the case where a; < < bi, it may be possible to split X% X • • • x X m . 
In this way, a more informative outcome may be obtained for one or both of the results 
of splitting. Such splits can be repeated as long as possible and necessary. 

3.6 Interval CSPs 

Problems in a wide variety of application areas can be expressed as CSPs. Domains can 
be as different as booleans, integers, finite symbolic domains, and reals. In this paper 
we consider Interval CSPs (ICSPs), which are CSPs where all domains are intervals 
and all constraints are primitive. 

ICSPs are important because they encapsulate what can be efficiently computed; 
they represent Layer 2 in the software architecture of Figure ^ The layer is distinct 
from Layer 3 because in Equation Q there typically occur atomic formulas that contain 
function symbols. 

To emphasize the role of ICSPs as a layer of software architecture, we view them 
as a virtual machine, with a function that is similar to those for Prolog or Java. Just as 
a program in Prolog or Java is translated to virtual machine instructions, a system such 
as Equation Q can be translated to an ICSP, as described in a later section. 

The instructions of the ICSP level are DROs, one for each constraint. As an example 
of such an ICSP virtual machine instruction, consider the DRO for product constraint. 
It reduces the box [a, b] x [c, d] x [e, /] to the box that has the projections 

V([a,b]n([ej]/[c,d})) 
<p([c,d]n([e,f]/[a,b])) 
ip([e,f]n([a,b}*[c,d])) 

Here ip is the function that yields the smallest interval containing its argument. 

Of particular interest is the effect of the DRO when all variables have [— oo, +oo] as 
domain. For each argument, the domain after application of the DRO is defined as the 
default domain of that argument. Typically, default domains are [— oo, +00] . Notable 
exceptions include the constraint sin(x,y) (defined as y = s'm(x)), where the default 
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domain of y is [—1, 1]. Another is sq(x,y) (defined as y = x 2 ), where the default 
domain of y is [0, oo]. 

A difference with other virtual machines is that a program for the ICSP virtual 
machine is an unordered collection of DROs. Programs for other virtual machines are 
ordered sequences of instructions. In those other virtual machines, the typical instruc- 
tion does not specify the successor instruction. By default this is taken to be the next 
one in textual order. Execution of the successor is implemented by incrementing the 
instruction counter by one. 

The simplicity of the instruction sequencing in conventional virtual (and actual) 
machines is misleading. Many instruction executions concern untypical instructions, 
where the next instruction is specified to be another than the default next instruction. 
Examples of such untypical instructions are branches (conditional or unconditional) and 
subroutine j umps . 

In the ICSP virtual machine, the DROs are the instructions, and they form an un- 
ordered set. Instead of an instruction counter specifying the next instruction, there is the 
active set of GPA containing the set of possible next instructions. Instead of an instruc- 
tion or a default rule determining the next instruction to be executed, GPA selects in an 
unspecified way which of the DROs in the active set to execute. In this way, programs 
can be declarative: instructions have only meaning in terms of what is to be computed. 
How it is done (instruction sequencing), is the exclusive task of the virtual machine. 

3.7 A canonical form for nonlinear numerical inequalities 

Equation Q may have multiple occurrences of variables in the same formula. As there 
are certain advantages in avoiding such occurrences, we rewrite without loss of gener- 
ality the system in Equation Q to the canonical form shown in Figure|3] 
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, 3/2 , • ■ 


• , Vn) 
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, 3/2 , • ■ 


■ , Vn) 
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, 3/2 


■ , Vn) 
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allEq , ui,2 , • . • , tti.m) 

allEq (v P) i , Up, 2 , • ■ • , v P ,n p ) 

Fig. 3. A system of non-linear inequalities without multiple occurrences of variables. 
Instead, the set {yi, . . . , y n } is partitioned into equivalence classes Vi, . . . , V v where 
Vj is a subset {w^i, . . . , Vj <nj } of {t/i, . . . , y n }, for j e {1, ... ,p}. An allEq constraint 
asserts that its arguments are equal. 

In Figure [3] the expressions for the functions gi, . . . , gj. have no multiple occur- 
rences of variables. As a result, they have variables y%, . . . ,y n instead of x\, . . . , x m , 
with to < n as in Equation Q. This canonical form is obtained by associating with 



8 



each of the variables Xj in Equation {0 an equivalence class of the variables in Fig- 
ure This is done by replacing in Equation Q each occurrence of a variable by a 
different element of the corresponding equivalence class. This is possible by ensuring 
that each equivalence class is as large as the largest number of multiple occurrences. 
The predicate allEq is true if and only if all its real-valued arguments are equal. 

An advantage of this translation is that evaluation in interval arithmetic of each 
expression gives the best possible result, namely the range of the function values. At 
the same time, the allEq constraint is easy to enforce by making all intervals of the 
variables in the constraint equal to their common intersection. This takes information 
into account from all k inequalities. If the system in its original form as in Equation^ 
with multiple occurrences, would be translated to a CSP, then only multiple occurrences 
in a single expression would be exploited at one time. 

In the coming sections and without loss of generality, we will only consider expres- 
sions without multiple occurrences of variables. 

3.8 Translating nonlinear inequalities to ICSPs 

ICSPs represent what we can solve. They consist of atomic formulas without function 
symbols that, moreover, have efficient DROs. Equation Q exemplifies what we want 
to solve: it consists of atomic formulas typically containing deeply nested terms. 

The tree form of a formula We regard a first-order atomic formula as a tree. The unique 
predicate symbol is the root. The terms that are the arguments of the formula are also 
trees and they are the subtrees of the root. If the term is a variable, then the tree only 
has a root, which is that variable. A term may also be a function symbol with one or 
more arguments, which are terms. In that case, the function symbol is the root with the 
argument terms as subtrees. 

In the tree form of a formula the leaves are variables. In addition, we label every 
node that is a function symbol with a unique variable. 

Any constants that may occur in the formula are replaced by unique variables. We 
ensure that the associated domains contain the constants and are as small as possible. 

Translating a formula to an ICSP The tree form of a formula thus labeled is readily 
translated to an ICSP. The translation has a set of constraints in which each element is 
obtained by translating an internal node of the tree. 

The root translates to p(xo, . . . , x n -i) where p is the predicate symbol that is the 
root and xq, . . . , % n ~\ are the variables labeling the children of the root. 

A non-root internal node of the form f(to, ■ ■ ■ , in-i) translates to F(xq, . . . , x n -i, y), 
where 

- y is the variable labeling the node 

- xq, . . . , x n -i are the variables labeling the child nodes 

- F is the relation defined by F(ao, . . . , a n _i, v) iff v = /(ao, . . . , a n -i) for all 
clq, . . . , a„-i,v. 
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3.9 Search 

Propagation may terminate with small intervals for all variables of interest. This is rare. 
More likely, propagation leaves a large box as containing all solutions, if any. To obtain 
more information about possibly existing solutions, it is necessary to split an ICSP into 
two ICSPs and to apply propagation to both. An ICSP S' is a result of splitting an ICSP 
S if S' has the same constraints as S and differs only in the domain for one variable, 
say, x. The domain for x in S' is the left or right half of the domain for x in S. 

A search strategy for an ICSP is a binary tree representing the result of successive 
splits. Search strategies can differ greatly in the effort required to carry them to comple- 
tion. The most obvious search strategy is the greedy strategy: the one that ensures that 
all intervals become small enough by choosing a widest domain as the one to be split. 

This is a plausible strategy in the case where the ICSP has a few point solutions. In 
general, the set of solutions is a continuum: a line segment, a piece of a surface, or a 
variety in a higher dimensional space that has positive volume. In such cases we prefer 
the search to result in a single box containing all solutions. Of course we also prefer 
such a box to be as small as possible. The greedy search strategy splits the continuum 
of solutions into an unmanageably large number of small boxes. It is not clear that the 
greedy strategy is preferable even in the case of a few well-isolated point solutions. In 
general we need a search strategy other than the greedy one. 

A more promising search strategy was first described in the absolve predicate of 
the BNR Prolog system 1 6 1 and by 1 3 1, where it is called box consistency. 

The box consistency search strategy selects a variable and a domain bound. Box 
consistency uses binary search to determine a boundary interval that can be shown to 
contain no solutions. This boundary interval can then be removed from the domain, thus 
shrinking the domain. This is repeated until a boundary interval with width less than a 
certain tolerance is found that cannot be shown to contain no solutions. When this is the 
case for both boundaries of all variables, the domains are said to be box consistent with 
respect to the tolerance used and with respect to the method for showing inconsistency. 
When this method is interval arithmetic, we obtain functional box consistency. When it 
is propagation, then it is called relational box consistency [ 14 1. 

All we need to know about search in this paper is that greedy search and box con- 
sistency are both search strategies and that both can be based on propagation. Box 
consistency is the more promising search strategy. Thus we need to compare interval 
arithmetic and propagation as ways of showing that a nonlinear inequality has no solu- 
tions in a given box. This we do in section|4] 

4 Propagation with selective initialization 

Suppose we have a term that can be evaluated in interval arithmetic. Let us compare 
the interval that is the result of such an evaluation with the effect of GPA on the ICSP 
obtained by translating the term as described in section l3~8l 

To make the comparison possible we define evaluation of a term in interval arith- 
metic. The definition follows the recursive structure of the term: a term is either a vari- 
able or it is a function symbol with terms as arguments. If the term is an interval, then the 
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result is that interval. If the argument is function / applied to arguments, then the result 
is / evaluated in interval arithmetic applied to the results of evaluating the arguments in 
interval arithmetic. This assumes that every function symbol denotes a function that is 
defined on reals as well as on intervals. The latter is called the interval extension of the 
former. For a full treatment of interval extensions, see 111 1 11 2181 . The following lemma 
appears substantially as Theorem 2.6 in |7|. 

Lemma 1. Let the a term that can be evaluated in interval arithmetic. Let the variables 
oftbexi,..., x n . Let y be the variable associated with the root of the tree form oft. 
Let S be the ICSP that results from translating t, where the domains of Xi,...,x n 
are X\, . . . , X n and where the domains of the internal variables are [— oo, +oo]. After 
termination of GPA applied to S, the domain for y is the interval that results from 
interval arithmetic evaluation oft. 

Proof. Suppose that a variable of a constraint has domain [— oo, +oo]. After applying 
the DRO for that constraint, this domain has become the result of the interval arithmetic 
operation that obtains the domain for this variable from the domains of the other vari- 
ables of the constraint. 

According to II 3111 . every fair sequence of DROs converges to the same domains for 
the variables. These are also the domains on termination of GPA. Let us consider a fair 
sequence that begins with a sequence s of DROs that mimics the evaluation of t in in- 
terval arithmetic. At the end of this, y has the value computed by interval arithmetic. 
This shows that GPA gives a result that is a subinterval of the result obtained by interval 
arithmetic. 

GPA terminates after activating the DROs in s. This is because in the interval arithmetic 
evaluation of t an operation is only performed when its arguments have been evaluated. 
This means that the corresponding DRO only changes one domain. This domain is the 
domain of a unique variable that occurs in only one constraint that is already in the 
active set. Therefore none of the DRO activations adds a constraint to the active set, 
which is empty after s. 

GPA yields the same result whatever the way constraints are selected in the active 
set. Therefore GPA always gives the result of interval arithmetic evaluation. However, 
GPA may obtain this result in an inefficient way by selecting constraints that have no 
effect. This suggests that the active set be structured in a way that reflects the structure 
of t. This approach has been taken in 141 141 . 

The proof shows that, if the active set had not contained any of the constraints only 
involving internal variables, these constraints would have been added to the active set 
by GPA. This is the main idea of selective initialization. By initializing and ordering 
the active set in a suitable way and leaving GPA otherwise unchanged, it will obtain the 
interval arithmetic result with no more operations than interval arithmetic. This assumes 
the optimization implied by the Totality Theorem in 1 10 1. 

Definition 1. A constraint is a seed constraint iff at least one of its variables has a 
domain that differs from the default domain assigned to that variable. 

For example, the term sin(a;i) + sin^) translates to an ICSP with constraints 
sum(u,v,y), sin(xi,u), and sin(x2,v). When the domains are [— oo. +oo] for all 



11 



variables, then the seed constraints are sin(xx, u) and sin(x2, v). When the domains 
are [—00, +00] for x\, X2, and y; [—1, 1] for u and v, then sum(u, v, y) is the one seed 
constraint. When the domains are [— 00, +00] for X2 and y; [— 1, 1] for x\, u and v, then 
the seed constraints are sum(u, v, y) and sin(xi, u). 

Definition 2. Let P SI (Propagation with Selective Initialization) be GPA except for the 
following modifications. 

(a) P SI only applies to ICSPs generated by translation from an atomic formula. 

(b) The active set is a priority queue that is ordered according to the distance from the 
root of the node that generated the constraint. The greater that distance, the earlier the 
item is removed from the queue. 

(c) The active set contains all seed constraints and no other ones. 

Lemma [J says that GPA simulates interval arithmetic as far as the result is con- 
cerned. It does not say anything about the efficiency with which the result is obtained. 
Theorem^says that PSI obtains the result as efficiently as it is done in interval arith- 
metic. This assumes the functionality optimization in the DROs 1 10 1 . 

Theorem 1. Let S be the ICSP obtained by translating a term t in variables x\, . . . , x n , 
where these variables have domains X\ , . . . , X n . Applying PSI to S terminates after 
selecting no constraint more than once. Moreover, the root variable ends up with Z as 
domain where Z is the interval resulting from evaluating t with x%, . . . ,x n substituted 
by Xi, . . . ,X n . 

Proof. Suppose GPA is applied to S in such a way that all non-seed constraints are se- 
lected first. The execution of the DRO corresponding to these non-seed constraints does 
not change any domains. Therefore these DRO executions do not add any constraints. 
As a result, the effect of applying GPA is the same as when the active set would have 
been initialized with only the seed constraints. 

Suppose the seed constraints are selected according to priority order. This ensures that 
no future constraint selection re-introduces a constraint previously selected. Thus GPA 
terminates after activating every seed constraint exactly once. 

Such an execution of GPA coincides step by step with that of PSI. As GPA terminates 
with the correct result, so does PSI. 

5 Using ICSPs to solve inequalities 

We briefly reviewed how interval arithmetic can solve systems of nonlinear inequalities. 
The fundamental capability turned out to be that of evaluating a term in interval arith- 
metic. We saw that this can also be done by applying propagation to ICSPs generated 
from arithmetic terms. We now investigate how to extend this to the use of ICSPs to 
solve nonlinear inequalities. 

5.1 Using ICSPs to solve a single inequality 

Suppose S is the ICSP resulting from translating 

gi(xx, . . . ,x m ) < 
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Let y be the variable labeling the left child of the root; that is, the variable representing 
the value of the left-hand side. Let X\, . , . , X m be the domains in S of xi, . . . , x m , 
respectively. 

Now suppose that GPA is applied to S. One possible initial sequence of DRO acti- 
vations is the equivalent of interval arithmetic evaluation of the left-hand side, leaving 
y < as the only constraint in the active set with the domain for y equal to [cjj, £>»], the 
value in interval arithmetic of gt(Xi, . . . , X m ). 

At this stage the DRO for y < is executed. If < dj, then failure occurs. If 
bi < 0, then the domain for y is unchanged. Therefore, no constraint is added to the 
active set. Termination occurs with nonfailure. There is no change in the domain of any 
of x%, . . . , x m . The third possibility is that a, < < bi. In this case, the domain for 
y shrinks: the upper bound decreases from bi to 0. This causes the constraints to be 
brought into the active set that correspond to nodes at the next lower level in the tree. 
This propagation may continue all the way down to the lowest level in the tree, resulting 
in shrinking of the domain of one or more of x\ , . . . , x m . 

Let us compare this behaviour with the use of interval arithmetic to solve the same 
inequality. In all three cases, GPA gives the same outcome as interval arithmetic: fail- 
ure or nonfailure. In the first two cases, GPA gives no more information than interval 
arithmetic. It also does no more work. 

In the third case, GPA may give more information than interval arithmetic: in addi- 
tion to the nonfailure outcome, it may shrink the domain of one or more of x\, ■ ■ ■ , x m . 
This is beyond the capabilities of interval arithmetic, which is restricted to transmit in- 
formation about arguments of a function to information about the value of the function. 
It cannot transmit information in the reverse direction. To achieve this extra capability, 
GPA needs to do more work than the equivalent of interval arithmetic evaluation. 

In the above, we have assumed that GPA magically avoids selecting constraints in 
a way that is not optimal. In such an execution of GPA we can recognize two phases: 
an initial phase that corresponds to evaluating the left-hand side in interval arithmetic, 
followed by a second phase that starts with the active set containing only the constraint 
y < 0. When we consider the nodes in the tree that correspond to the constraints that 
are selected, then it is natural to call the first phase bottom-up (it starts at the leaves 
and ends at the root) and the second phase top-down (it starts at the root and may go 
down as far as to touch some of the leaves). The bottom-up phase can be performed 
automatically by the PSI algorithm. 

The start of the top-down phase is similar to the situation that occurs in search. In 
both search and in the top-down phase a different form of selective initialization can be 
used, shown in the next section. 

The bottom-up phase and the top-down phase are separated by a state in which the 
active set only contains y < 0. For reasons that become apparent in the next section, 
we prefer a separate treatment of this constraint: not to add it to the active set and to 
execute the shrinking of the domain for y as an extraneous event. This is then a special 
case of termination of GPA, or its equivalent PSI, followed by the extraneous event of 
shrinking one domain. 

The Pseudo-code for PSI algorithm is given in Figure|4] 

The correctness of PSI algorithm can be easily deduced from the following theorem. 



13 



let the active set A be a priority queue in which the constraints are 
ordered according to the level they occupy in the tree, 

with those that are further away from the root placed nearer to the front of the queue 
put only seed constraints into A 
while ( A ± 0) { 

choose a constraint C from A 

apply the DRO associated with C 

if one of the domains is empty, then stop 

add to A all constraints involving variables whose domains have changed, if any 
remove C from A 

} 

Fig. 4. Pseudo-code for Propagation with Selective Initialization (PSI). 

Theorem 2. Consider the ICSP S obtained from the tree T of the atomic formula 
gi(xi, . . . ,x n ) < 0. Suppose we modify GPA so that the active set is initialized to 
contain instead of all constraints only seed constraints. Suppose also that the active 
set is a priority queue in which the constraints are ordered according to the level they 
occupy in the tree T, with those that are further away from the root placed nearer to 
the front of the queue. Then GPA terminates with the same result as when the active set 
would have been initialized to contain all constraints. 

Proof. As we did before, suppose that in GPA the active set A is initialized with all 
constraints such that the seed constraints are at the end of the active set. Applying any 
DRO of a constraint that is not a seed constraint will not affect any domain. Thus, 
the constraints that are not seed constraints can be removed from the active set without 
changing the result of GPA. Since the GPA does not specify any order, A can be ordered 
as desired. Here we choose to order it in such a way we get an efficient GPA when used 
to evaluate an expression (see previous section). 

6 Selective Initialization for search 

Often we find that after applying GPA to an ICSP S, the domain X for one of the 
variables, say x, is too wide. Search is then necessary. This can take the form of splitting 
S on the domain for x. The results of such a split are two separate ICSPs S\ and S2 that 
are the same as S except for the domain of x. In Si, x has as domain the left half of X; 
in S2, it is the right half of X. 

However, applying GPA to Si and S2 entails duplicating work already done when 
GPA was applied to S. When splitting on x after termination of the application of 
GPA to S, we have the same situation as at the beginning of the downward phase of 
applying GPA to an inequality: the active set is empty and an extraneous event changes 
the domain of one variable to a proper subset. 

The following theorem justifies a form of the PSI algorithm where the active set is 
initialized with what is typically a small subset of all constraints. 
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Theorem 3. Let T be the tree obtained from the atomic formula gi(x±, . . . , x n ) < 0. 
Let S be the ICSP obtained from T. Let x be a variable in S. Suppose we apply GPA 
to S. After the termination of GPA, suppose the domain of x is changed to an interval 
that is a proper subset of it. If we apply GPA to S with an active set initialized with 
the constraints only involving x, then GPA terminates with the same result as when the 
active set would have been initialized to contain all constraints. 

Proof. To prove Theorem|3j we should show that initializing GPA with all constraints 
gives the same results as when it is initialized with only the constraints involving x. 

Since no ordering is specified for the active set of GPA, we choose an order in 
which the constraints involving x are at the end of the active set. Because DROs are 
idempotent, all constraints at the front of the active set, different from those involving 
x, do not affect any domain. Thus removing them from the active set in the initialization 
process does not change the fixpoint of the GPA. Thus, Theorem[3]is proved. 

7 Further work 

We have only considered the application of selective initialization to solve a single 
inequality. A conjunction of inequalities such as Equation Q can be solved by solving 
each in turn. This has to be iterated because the solving of another inequality affects 
the domain of an inequality already solved. This suggests performing the solving of 
all inequalities in parallel. Doing so avoids the waste of completing an iteration on the 
basis of unnecessarily wide intervals. It also promises speed-up because many of the 
DRO activations only involve variables that are unique to the inequality. In the current 
version of the design of our algorithm, we combine this parallelization with a method 
of minimizing the complexity usually caused by multiple occurrences of variables. 

8 Conclusions 

Before interval methods it was not clear how to tackle numerically realistic optimization 
models. Only with the advent of interval arithmetic in the 1960s 1 1 1 1 one could for the 
first time at least say: "If only we had so much memory and so much time, then we 
could solve this problem." 

Interval arithmetic has been slow in developing. Since the 1980s constraint pro- 
gramming has added fresh impetus to interval methods. Conjunctions of nonlinear in- 
equalities, the basis for optimization, can be solved both with interval arithmetic and 
with constraint programming. In this paper we relate these two approaches. 

It was known that constraint propagation subsumes interval arithmetic. It was also 
clear that using propagation for the special case of interval arithmetic evaluation is 
wasteful. In this paper we present an algorithm for propagation by Selective Initializa- 
tion that ensures that propagation is as efficient in the special case of interval arithmetic 
evaluation. We also apply Selective Initialization for search and for solving inequali- 
ties. Preliminary results on a parallel version of the methods presented here suggest that 
realistic optimization models will soon be within reach of modest computing resources. 
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